Ab initio and finite-temperature molecular dynamics studies of lattice 

resistance in tantalum 
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This manuscript explores the apparent discrepancy between experimental data and theoreti- 
cal calculations of the lattice resistance of bcc tantalum. We present the first results for the 
temperature dependence of the Peierls stress in this system and the first ab initio calculation 
of the zero-temperature Peierls stress to employ periodic boundary conditions, which are those 
best suited to the study of metallic systems at the electron-structure level. Our ab initio value 
for the Peierls stress is over five times larger than current extrapolations of experimental lat- 
tice resistance to zero-temperature. Although we do find that the common techniques for such 
extrapolation indeed tend to underestimate the zero-temperature limit, the amount of the un- 
derestimation which we observe is only 10-20%, leaving open the possibility that mechanisms 
other than the simple Peierls stress are important in controlling the process of low temperature 
slip. 



I INTRODUCTION 

The study of the plasticity of crystalline materials is 
a rich many-body problem involving physics on multiple 
length scales, with many remaining unexplained myster- 
ies. The plasticity of bcc metals, for instance, is particu- 
larly challenging. Unlike their fee and hep counterparts, 
the bcc metals exhibit many active slip planes, have a 
strong temperature dependence in their plasticity and 
violate the simple empirical Schmid law.^ Moreover, the- 
oretical calculations of the most basic question in plastic- 
ity, the stress needed to induce yield at low temperature 
in a pure sample, differ from experimental extrapolations 
by over a factor of two.^ The purpose of this work is to 
provide needed insight into this discrepancy. 

Ultimately, it is the physics of the (111) screw disloca- 
tion defect which controls the low-temperature plasticity 
of bcc materials.^' ^ The Peierls stress, the yield stress 
at which these dislocations first begin to move sponta- 
neously, is difficult to compare directly with experiment. 
Whereas most computational work on the Peierls stress 
measures the stress to move an isolated, infinitely straight 



dislocation at zero temperature 
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experiments mea- 



sure the Peierls stress at a finite temperature in systems 
with many interacting, curved dislocations and in media 
with defects and surfaces. As an example of the present 
challenges which the literature faces, using model gener- 
alized pseudopotential theory (MGPT),^^ Yang et. al}^ 
predict for the T=0 Peierls stress a value 2.5 times greater 
than experimental extrapolations.*^ Because such poten- 



tials are not based upon first principles, it is impossible 
to determine a priori whether this discrepancy is due to 
the inter-atomic potential, the environmental complexi- 
ties discussed above, or to a flaw in our understanding 
of the relation between the Peierls stress and the experi- 
ments. 

As it is a daunting experimental task*^ to observe 
properties of a single dislocation locked deep in the heart 
of a material, accurate theoretical calculations of such 
systems is essential. Nearly all theoretical calculations to 
date, concerning such dislocations, have relied upon em- 
pirical potentials. *°"*^' Given the empirical na- 
ture of such calculations, the complex directional bonding 
properties of bcc materials, and the lack of direct com- 
parison with experiments for validation, first principles 
ab initio calculations of dislocations in such systems are 
clearly needed. Ismail-Beigi and Arias^° were the first 
to show that density functional theory calculations were 
crucial in understanding the fundamental properties of 
the (111) screw dislocation core structure in bcc molyb- 
denum and tantalum. Until that work, most computa- 
tional studies based on empirical potentials, ii^i^ gyp_ 
ported the idea that the dislocation core breaks symme- 
try, with two energetically equivalent ground state struc- 
tures which spread outward along two different equiva- 
lent sets of three {110} planes,'^ similar to the concept 
originally proposed by Hirsch.^*'^'^ Until the availability 
of the ab initio calculations, the prevailing view of the 
violation of the Schmid law in the bcc metals was based 
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upon this structure.^ Ismail-Beigi and Arias, ^'^ in con- 
trast, showed that for both molybdenum and tantalum 
the ground-state structure within density functional the- 
ory was a non-degenerate symmetric core, strongly sup- 
porting the work of Suzuki and Takeuchi^^'^^'^'' which 
first suggested that it is the Peierls potential itself that 
controls the lattice resistance and not the details of the 
core structure. To help resolve the discrepancy between 
theoretical and experimental Peierls stresses, the work 
below provides a reliable ab initio prediction of the Peierls 
stress in bcc tantalum which is free of the unrealistic elec- 
tronic boundary conditions employed in the only other ab 
initio prediction of the Peierls stress.^' ^ 

Here, we show that the Peierls stress, calculated within 
density function theory, is over a factor of five larger than 
expected from extrapolation of experimental results.^® 
This supports the view that the discrepancy between the 
experimental and computational predictions are largely 
due to the aforementioned environmental complexities, 
to a flaw in relating the experimental data to the Peierls 
stress, or to a combination of both. 

To further explore possible physical effects leading to 
this discrepancy, we study the extrapolation of experi- 
mental data to determine the zero-temperature Peierls 
stress. Such extrapolations generally employ fits from 
mesoscopic or thermodynamics/kinetic models.^''' 
However, it has not been established that such models 
can accurately describe the lowest temperature regime 
correctly, placing doubt on the quality of these extrapola- 
tions. To address this issue, the work below also provides 
the first temperature- and orientation- dependent study 
of the Peierls stress in a bcc metal. We moreover show 
that extrapolation of our finite temperature results using 
a current fitting model leads to an underestimation of the 
zero-temperature Peierls stress. This underscores the dif- 
ficulty in extrapolating the experimental data accurately 
but does not fully account for the observed discrepancy. 

In Section II, this manuscript reviews the various tech- 
niques in use for calculation of the Peierls stress in the 
context of efficacy for application to ab initio calcula- 
tions. Section III gives the first calculation of the tem- 
perature and orientation dependent Peierls stress in a bcc 
material. Section IV describes our technique for obtain- 
ing Peierls stresses within small unit cells with periodic 
boundary conditions. Finally, Section V presents our ab 
initio prediction for the Peierls stress and compares and 
contrasts it to currently available experimental and com- 
putational values. 



II BOUNDARY CONDITIONS 

The fundamental distinction among theoretical ap- 
proaches to calculation of the Peierls stress is the choice 
of boimdary condition. The literature describes three 
types of boundary conditions: cylindrical boundary 
conditions,^"^' Greens function (or "flexible") 



boundary conditions,"^' ^^^"'^^'^^ and periodic boundary 
conditions.^'®' We now briefly review each with em- 
phasis on the unique challenges of ab initio electronic 
structure calculations. 



A Cylindrical Boundary Conditions 



In the practice of cylindrical boundary conditions, 

anisotropic elasticity theory"^""''^ is used to generate a dis- 
location in the center of a cylinder. The cylinder is then 
separated into inner and outer regions. The atoms in the 
outer region are held flxed to the solution of anisotropic 
elasticity theory while the atoms in the inner region re- 
lax under the inter-atomic forces. To calculate the Peierls 
stress, a stress is applied to the system until the disloca- 
tion moves. 

This approach suffers numerous draw backs when ap- 
plied to density functional theory. To avoid surface ef- 
fects and to properly account for the non-linear nature 
of the dislocation, such cylinders generally have to be 
quite large. First, even the outer cylinder is of finite size 
and therefore the outer region must be sufficiently large 
enough so that forces generated by it onto the inner re- 
gion are equivalent to those generated from an infinite 
continuum. The inner region also must be sufficiently 
large to mitigate two effects. The inner region iramt be 
large enough so that linear elasticity theory represents 
well the forces which it imposes on the outer region. The 
inner region also must be large enough so that motion 
of the dislocation is not adversely affected by the fixed 
outer region, which is a concern because the fixed outer 
region represents the displacement field when a disloca- 
tion is at its center and therefore generates a extraneous 
force which tends to prevent motion of the dislocation.^^ 
When using simple, inter-atomic potentials, the use of 
large cylinders mitigates all of these effects. However, 
this approach is not viable for density functional calcu- 
lations with their extreme computational demands. 

This approach, moreover, is particularly ill-suited for 
electronic structure calculations because the artificial 
surface at the outside of the outer region, being far dif- 
ferent from the bulk, give rise to strong scattering of 
the electrons far different than would an infinite contin- 
uum. This is particularly problematic for metals, be- 
cause the single-particle density matrix, which quantifies 
the effects of this scattering on the inter-atomic forces, 
decays only algebraically in metals. The following sub- 
section demonstrates that the boundary regions should 
be quite large in order for these surface effects not to re- 
sult in large fictitious forces in the active region of the 
calculation. 
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B Greens Function Boundary Conditions 

The use of Greens function, or flexible, boundary con- 
ditions,^' ^^'^^ is an effective way to reduce the size of the 
simulation cell. This approach also employs a cylindrical 
geometry. However, rather than the "inner" and "outer" 
atomistic regions of the cylindrical boundary approach, 
the Greens function approach employs three inter-atomic 
regions: an inner "core" region containing the center 
of the dislocation, an intermediate "buffer" region, and 
an outermost "continuum-response" region. With proper 
implementation, the outer and inner regions couple only 
indirectly through the response of the buffer region. 

In this method, all three regions respond to the pres- 
ence of a dislocation; however the response of each region 
is treated differently through a number of steps. Initially, 
all regions are displaced by the solution to anisotropic 
elasticity theory. Each iteration then begins by relaxing 
the atoms in the core region according to the forces which 
they experience, as computed from cither an inter-atomic 
potential or an ab initio method. The forces generated 
from the mismatch between the outer and inner regions, 
which the cylindrical approach above ignores, are then 
relieved by moving the atoms of all three regions accord- 
ing to the elastic Greens function solution, leaving only 
the nonlinear effects from the core region unaccounted. 
The next iteration then begins by relaxing these forces as 
described above. Iterations proceed until until the forces 
in the core and buffer region are negligible. 

What distinguishes this approach from simple cylin- 
drical boundary conditions is that the continuum region, 
via the Greens function response, is allowed to respond 
to the motion of the dislocation and to the elastic re- 
sponse generated by the core region as the dislocation 
moves. So long as the continuum response region (a) ac- 
curately represents the structure induced by the presence 
of the dislocation and (b) is sufficiently wide to properly 
reproduce the forces on the atoms in the buffer and inner 
regions, this approach accurately describe basic proper- 
ties of a dislocation. 

In order for the first assumption (a) above to hold, the 
inner core region must be sufficiently large to contain all 
atoms with displacements outside of the linear regime 
and the buffer region must be sufficiently wide so that 
displaced atoms in the core have no effect on the forces 
experienced in the continuum-response region. The sec- 
ond assumption (b) requires that the continuum-response 
region to be sufficiently large so that its termination has 
no effect on the forces on the atoms in the buffer or inner 
region. The radius of the calculation must therefore ex- 
ceed the sum of the non-linear core radius plus twice the 
range over which motion of atoms creates forces within 
the lattice. As the latter range can be quite large for 
electronic structure calculations in metals, the applica- 
tion of this approach to electronic structure calculations 
can be problematic. 

The Greens function approach has predicted success- 



fully dislocation properties when applied to time con- 
suming empirical potentials^' which have a limited in- 
teraction range. The approach also has been applied to 
density functional calculations of the Peierls stress for 
molybdenum and tantalum,^' ^ where its application is 
more questionable due to the above interactions. In these 
latter works, the artificial boundary on the outside of the 
continuum-response region have been treated in one of 
two ways,**'*' either by keeping the surface free in vac- 
uum or by embedding in periodic boundary conditions 
with the vacuum filled with material which must contain 
severe domain boundaries due to the incompatibility of 
a net Burgers vector with periodic boundary conditions. 

To gauge the effects which this artificial boundary 
may have and how far these effects penetrate from the 
continuum-response region into the buffer region, we per- 
form a test calculation within the density functional the- 
ory pseudopotential approach"^^ of the magnitude of the 
forces generated onto the system due to the presence of 
a domain boundary similar to those in the works cited 
above.®' ^ 

For this calculation, we employ the same computa- 
tional procedure as for our production calculations in 
Section V. Here, however, as this is a test, we employ 
only a single fc-point to sample the Brillouin zone (F). 
We begin with an orthorhombic cell of 24 atoms of tan- 
talum in a bulk arrangement with supercell lattice vec- 
tors fi = a[110], f2 = 4a[112] and = a/2[lll]. We 
choose this cell because its length along r2 is the same as 
the smallest simulation cell used in References 8 and 9. 
We then generate a domain boundary at the edge of cell 
along the (112) plane by changing the lattice vector f2 to 
r2 = 4a[112] -l-ars and holding the atoms in the unit cell 
fixed in their bulk locations, a is chosen such that the 
shift is small and the nearest neighbor distance is always 
within 95% of the bulk, representing even less of a distur- 
bance that in Reference 8, where atom were within 90% 
of the bulk nearest neighbor distance. To estimate the 
effect of the scattering of electrons at the domain bound- 
ary on the inter- atomic forces, we hold the atoms fixed 
and compute the ab initio forces acting upon them. 

Figure 1 shows the forces along the [111] direction as 
a function of distance from the center of each domain. 
Note that relatively large forces develop deep within the 
cell. This data indicates that the continuum-response re- 
gion should be quite large 5 — 10 A) in order for the 
response of the electrons not to adversely effect the forces 
in the buffer region. Note also that the the buffer region 
should be of similar width to prevent forces from the non- 
linear displacements in the core from penetrating into the 
linear continuum-response region. Such large continuum- 
response and buffer regions can make the calculation in- 
feasible with current computational techniques. 

In fact, the only density functional calculations of the 
Peierls stress in this system to date employ the Greens 
function method but with a distance from the buffer re- 
gion to the domain boundary of only ^ 3.7 A. It thus 
is unclear whether the continuum region in these calcu- 
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lations is sufficiently large to lead to reliable results and 
clearly further calculations are needed to support those 
results. Below, we provide just such calculations using 
the method of periodic boundary conditions, which per- 
turb the electronic system far less than the introduction 
of domain boundaries. 
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FIG. 1: Force on the atoms, along the [111], due to the 
presence of a domain boundary. The forces are plotted 
as a function of distance from the center of the unit cell. 
The domain boundary is generated such that the nearest 
neighbor distance is always within 95% of the bulk (a = 
1/4). 



C Periodic Boundsiry Conditions 

The final common choice for boundary conditions is 
to repeat the dislocation core periodically throughout 
space, so that the dislocation is no longer isolated, but 
embedded in bulk material containing an array of dislo- 
cations. Consistency with periodic boundary conditions 
then demands that the unit cell contain an even number 
of dislocations with Burgers vectors of alternating sign, 
arranged typically in dipolar'^^"'^ or quadrupolar^°' ^ 
arrays. For static properties of the dislocation core, such 
cells give reliable results as the elastic fields of the sur- 
rounding dislocations effectively cancel at the location of 
each core. 

Periodic boundary conditions can also been used to 
calculate the Peierls stress.^' Care must be taken, 
however, because the symmetry of the dipole/quadrupolc 
array breaks as the dislocations move and the dislocation- 
dislocation interaction then ceases to be negligible. The 
use of large unit cells can control this effect;^'^ however, 
such a direct, brute-force approach is not practical for 
computationally demanding ab initio calculations. To 
make such calculations feasible, they must occur in small 
periodic cells, thereby demanding proper accounting of 
the dislocation-dislocation interactions. 



We have shown in another work^ that, under certain 
conditions, such interactions can be accounted accurately 
with minimal extra computational effort, so that accurate 
values of the Peierls stress can be obtained from density 
functional calculations in periodic cells. As the residual 
errors with this approach are associated with the bound- 
ary conditions, the magnitude of such an error can be 
tested by other computational methods, such as empir- 
ical potentials. Our previous work shows that that this 
residual error is relatively small, ^ a fact which we confirm 
explicitly below. 

Because the deviations from the bulk arrangement 
at periodic boundaries are relatively mild, such calcu- 
lations are ideal for mitigating electronic boundary ef- 
fects. Given the simplicity of working with these bound- 
ary conditions and the possibility of the extraction of ac- 
curate values for the Peierls stress from small unit cells, 
we choose to work with periodic boundary conditions 
throughout this work. Section IV outlines our procedure 
for calculating the Peierls stress while working with peri- 
odic boundary conditions and describes the sources and 
the magnitude of the residual errors. (See Reference 5 
for a full discussion of these issues.) 



Ill DEPENDENCE OF THE PEIERLS STRESS 
ON ORIENTATION AND TEMPERATURE 

To illustrate the complexities of relating computational 

predictions to experimental findings, we now explore the 
dependence of the Peierls stress in bee tantalum on orien- 
tation and temperature. The strong dependencies which 
we shall find underscore the unique properties of dislo- 
cations in bee metals. To clarify, as some authors use 
slightly different definitions for the Peierls stress, here 
we consider the Peierls stress as the value of the stress 
on the maximum resolved shear stress plane (maximum 
value of the shear stress along the [lll]-direction) when 
the dislocation first moves to a different equilibrium po- 
sition. 

Despite recent advances in ab initio quantum me- 
chanical methods, such methods are still too computa- 
tionally intensive to study such properties as the tem- 
perature dependence of the Peierls stress. Therefore, 
for these calculations, we employ a molecular dynam- 
ics (MD) framework carried out using a first-principles- 
based, many body force field (FF) for tantalum, which 
we denote qEAM, which we have developed to allow ac- 
curate and computationally eflacient evaluation of atomic 
interactions.^' '^^ 

As described above, we carry out these calculations 
within periodic boundary conditions. The super-cell con- 
sists of a quadrupolar arrangement"^''''^* of dislocations 
containing 5670 atoms with lattice parameters Ox = 
70.59 A, Oy = 73.39 A and = 20.11 A, where the 
X-, y-, and z- axes of our coordinate system are along 
[110], [112], and [111] directions, respectively. As we have 
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shown in another work,^ such a cell gives very accurate 
values for the Peierls stress. 



A Orientation dependence of the zero-temperature 
Peierls stress of the (111) screw dislocation 

To calculate the zero temperature Peierls stress, we 
start with a fully relaxed quadrupole dislocation configu- 
ration at zero stress and increase the stress in steps of 50 
MPa until the dislocations move. Once the dislocations 
move, we restart the calculation from the structure equi- 
librated just prior to the motion and increase the stress 
in smaller steps (5 MPa) in order to more narrowly define 
the critical stress. At each incremental target stress, we 
relax the atoms and stresses in the cell by running two 
very low temperature (T=0.001 K) MD simulations. The 
first run is for 15 ps at constant stress and temperature 
(NctT ensemble) using a Rahman-Parrinello barostat*° 
and a Hoover'*^ thermostat, and the second run is for 50 
ps at constant volume and temperature (NVT ensemble). 
We find this approach to be quite stable for relaxing the 
cell and the atoms of the system. 

The (111) screw dislocation has three equivalent {112} 
and three equivalent {110} potential slip planes, with 
such planes occurring at 30° intervals. To study the ori- 
entation dependence of the Peierls stress, we apply three 
types of pure shear stress to the system: a a^z stress, 
a positive ayz stress and a negative stress. (Note 
that with coordinate axes as defined above, the z-axis lies 
along the dislocation line.) These stresses leads to forces 
on the dislocation in the (112), (llO)-twinning and (110)- 
antitwinning directions, respectively.^'^ Along these di- 
rections, we find Peierls stresses of rii2 = 655 MPa, 
Ttwin = 575 MPa, and Tantitwin = 1075 MPa, respectively. 
Table I shows that our essentially zero-temperature re- 
sults arc in good agreement with those of Yang and col- 
laborators,^*^ who employed model generalized pseudopo- 
tential theory (MGPT),^'^ a different inter-atomic poten- 
tial. The result of such a strong dependence of the Peierls 
stress on orientation is consistent with the experimentally 
observed breakdown of the Schmid law in bcc metals. 



Potential 


Twin 


(112) 


Anti-Twin 


Asymmetry 


qEAM FF 

mgpt"' 


575 MPa 
605 MPa 


655 MPa 
640 MPa 


1075 MPa 
1400 MPa 


1.6412 
2.29 



TABLE I: Peierls stress for the (111) screw dislocation 
in tantalum in the twinning, (112) and anti-twinning di- 
rections; the last column shows the ratio between anti- 
twinning and twinning Peierls stresses. MGPT results 
from Yang et al.^'' 



To make quantitative comparison with experiments, 
which are carried out at nonzero temperature, we com- 
pare our results to those of Tang et. al..^^ who fit ex- 
perimental data**^ to a mesoscopic model and then ex- 
trapolate to extract the zero-temperature Peierls stress. 



Their predicted value of 248 MPa for the (112) Peierls 
stress is over a factor of two lower than our result. This 
type of discrepancy, where the theoretical Peierls stress 
overestimates the zero-temperature extrapolation of the 
experimental data by a factor of two to three, is quite gen- 
erally observed.^ This discrepancy may be due either to 
inaccuracies in the theoretical calculations or, perhaps, to 
a flaw in the comparison between the zero-temperature 
extrapolation of the experimental data and theoretical 
predictions. 



B Temperature dependence of the Peierls stress of 
the (111) screw dislocation 

To explore potential difficulties with the zero- 
temperature extrapolation, we now present what to our 
knowledge is the first temperature-dependent study of 
Peierls stress using a realistic potential for a bcc metal. 
For these calculations, wc continue to employ the qEAM 
FF and begin with the zero-temperature, equilibrated 
structures. We then apply various constant shear stresses 
(lower than the T=0 Peierls stress) to the system while 
slowly increasing the temperature (in steps of 10 K) until 
the dislocations move. Similarly to the T = 0.00 IK case, 
for each temperature we first run for 10 ps in the NuT 
ensemble and then for 25 ps in NVT ensemble. 

Because the Peierls stress can depend on the rate at 
which the strain is applied, to place our results in con- 
text, we first estimate the strain rate in our computa- 
tions. The strain rate is approximately 7 = pvdh, where 

is the dislocation velocity, p is the dislocation density 
and h is the Burgers vector. Using a dislocation den- 
sity typical of the experiments^® (p = 10^^1/m^) and 
estimating the dislocation velocity as the ratio between 
the distance traveled in one jump (l/3o(112) = 2.717 A) 
and the simulation time (35 ps), we obtain an effective 
strain rate of ^ 10^1/s, which is large compared to the 
strain rates (4 x 10~^) in the experiments used for the 
zero-temperature extrapolations.^^' 

Figure 2 summarizes our results for the temperature 
dependence of the Peierls stress as a function of tem- 
perature for the three directions ((112), twinning and 
anti-twinning). As expected, the Peierls stress obtained 
from our MD simulations decreases rapidly with increas- 
ing temperature, particularly for very low temperatures. 
It is important to mention that, although our simulations 
are three dimensional, the dislocations move as straight 
lines without the formation of double kinks because our 
simulation cell is only seven Burger's vectors long along 
the dislocation lines. Such double kinks are quite impor- 
tant at finite temperatures as they tend lower the lattice 
resistance at non-zero temperatures. Our results arc ap- 
proximately a factor of 2 — 4 larger than the fit of Tang 
et. al.^^ to the experimental data of Wasserbach.''^ We 
feel that this is reasonable, considering the facts that our 
simulation cells do not allow for double kink formation 
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and that, as discussed above, our strain rates are much 
higher than those in the experiments. ^^'^^ 

Experimental extrapolations of the Peierls stress to 
zero temperature generally come from mesoscopic or ki- 
netic/thermodynamic model^®'^^"^'' fits to experimental 
data and extrapolated to zero temperature. To explore 
the effects of this procedure, we fit our atomistic data to 
such a model, perform the extrapolation and then com- 
pare with our direct zero-temperature results. 



1400 I ■ , , ■ r- 

1200 - Ta qEAM FF 

1/2a<1 11 > screw dislocafon 



„ 1000 i 

Q_ 




I ■ ' ■ ' ■ ■ ■ ■ ■ ' ■ 1 

50 100 150 200 250 300 

Temperature (K) 



two consecutive Peierls valleys. Physically, E is the 
minimum energy to form a double kink, L^"^^ is the min- 
imum length for this double kink and Tq is the stress, 
whose work to move a dislocation a distance Ip is equal 
to E^'""^. 

Figure 2 shows the fit of Equation (1) to our atomistic 
data. To mimic how zero-temperature lattice resistances 
are generally extracted, we adjust the three unknown pa- 
rameters (to, E^^'^^ and io'"'') to fit our higher tempera- 
ture data only. Intriguingly, extrapolation of our higher 
temperature results to zero-temperature leads to an un- 
derestimation of the Peierls stress of between 10% and 
20%. We would also expect that a fit to data from a 
cell sufficiently large to allow for double kink formation 
(which is active in the experiments at non-zero temper- 
atures) would lead to an even a larger underestimation 
of the zero-temperature stress. These results therefore 
suggest that the general discrepancy between extrapo- 
lated experimental values and the calculated values for 
the zero-temperature Peierls stress may be the result of 
failure of nonzero-temperature models to describe prop- 
erly the low-temperature regime. This illustrates one 
possible difficulty in relating the experiments to compu- 
tational predictions and underscores the need for first 
principles studies of such a phenomena. 



FIG. 2: Temperature dependence of the Peierls stress 
along various directions: the (112), twinning and anti- 
twinning directions. The fits are done for the high tem- 
perature data using Equation (1). The temperature is in 
Kelvin, and the stress in MegaPascals. 

For the fit we use an analytical expression for the de- 
pendence of the Peierls stress with temperature at con- 
stant strain rate^^' based on a mechanism involving 
double kink nucleation and propagation. This model 
gives for the temperature dependent Peierls stress 

where fi is l/ksT, with ks being Boltzmann's constant 
and T the absolute temperature, E^™"^ is the kink en- 
gj.gy 30 jg ^YiQ effective Peierls stress and 7o™'' 
reference strain rate. Here, the effective Peierls stress is 

^kink 

^0 = ^Zi^ 
and the reference strain rate is 

70^-'^ = 2hplpvn, (3) 

where h is the Burgers vector, L^™^ is the kink length, p 

is the dislocation density and is the attempt frequency 
which may be identified with the Debye frequency to a 
first approximation,^^' and Ip is the distance between 



IV ACCURATE PEIERLS STRESS 
CALCULATIONS IN SMALL PERIODIC 
CELLS 

Having underscored the need for first principles elec- 
tronic structure studies and already determined the most 
effective boundary conditions for such studies as periodic, 
we now focus on determination of the cell of minimal size 
appropriate to calculation of the zero-temperature Peierls 
stress for a (111) screw dislocation in a bcc metal when 
the maximum resolved shear stress is along a {110} plane. 

To minimize image effects, we employ periodic bound- 
ary conditions with a quadrupolar unit cell. Figure 3 
shows a differential displacement map'^ of such a cell of 
size 42 Ax 41 A in the plane perpendicular to the Burgers 
vector. In such a map, the dots indicate columns of atoms 
along the [111], and the vectors between the columns of 
atoms indicate the relative shift along the Burgers vec- 
tor due to the presence of a dislocation between each 
pair of columns, with the vectors scaled so that a vector 
of full length between the columns corresponds to 1/3 a 
Burgers vector. In this ground state structure of the dis- 
location, triads of full length vectors surrounds the center 
of each dislocation, corresponding to a new displacement 
by a full Burgers vector upon completion of a closed loop 
about the each center. In practice, because of symmetry, 
the quadrupole cell may be reduced in half containing two 
dislocations, when lattice vectors are properly chosen. 

As Section II notes, for calculations of static proper- 
ties such as the ground-state dislocation core structure, 
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the strain fields from the surrounding dislocations in a 
quadrupolar array essentially cancel at each dislocation 
core. However in a dynamical problem such as the calcu- 
lation of the Peierls stress, the dislocations begin to in- 
teract with the stress fields of the others as they begin to 
move. Our previous work^ shows that, for the particular 
geometry considered here, accurate values for the Peierls 
stress can be extracted from quite small unit cells pro- 
vided the proper procedure is followed. We now outline 
that procedure while reviewing the relevant background. 
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FIG. 3: Quadrupolar unit cell of size 42 Ax 41 A: relative 
displacements of neighboring columns of atoms (arrows), 
dislocation centers (plus signs). 



A Calculation of Peierls stress within periodic 
boundary conditions 

To calculate the Peierls stress in a small periodic cell, 
we begin with lattice vectors appropriate to bulk material 
in the absence of dislocations and then, while relaxing the 
internal coordinates of the cell, apply increasing pure e^z 
strains (same Cartesian coordinates from Section II C) 
until the dislocations move. Such a strain drives the 
dislocation along the [112] direction. Before the dis- 
location moves, the strain energy of the cell increases 
quadratically, 

E = lc'el,, (4) 

where C is an elastic constant associated with the 

quadrupole unit cell which can be extracted simply from 
the energies of the cell as the strain increases. The stress 
associated with this strain is 

(^xz = Ccxzi (5) 

SO that at the strain at which the dislocation moves. 
Equation (5) gives the Peierls stress. 



The great benefit of the above procedure is that it re- 
quires a minimal search through phase space in order 
to calculate the Peierls stress and accounts accurately 
for the effects of the dislocation-dislocation interactions. 
The C elastic constant, which comes as a direct byprod- 
uct of the procedure and requires no further calculations, 
suffices to account accurately for the leading order effects 
of the inter-dislocation interactions. This correction cap- 
tures most of the effects from working with small unit 
cells and can differ from that of an equivalent bulk cell 
by a value of two or more. 

The above procedure involves several approximations 
requiring justification. First, we apply strain relative to 
the lattice vectors of the bulk cell in the absence of the 
dislocation, rather than those of the quadrupole array. 
Previously, we have shown through explicit calculation 
on model inter-atomic potentials that working with the 
relaxed lattice vectors of the quadrupolar array docs not 
improve the value calculated for the Peierls stress and 
therefore is not needed.^ The reason for this is that al- 
though working with the bulk lattice vectors generates 
artificial stresses, these arc primarily diagonal (a^i), be- 
cause the greatest effect of the presence of the disloca- 
tions is to dilate the system.^ In the present geometry, 
such diagonal stresses result only in a constant shift in the 
energy over the range of applied strain and do not gener- 
ate driving (Pcach-Kohlcr'^°) forces on the dislocations. 

Some slight care must be taken with the above ar- 
gument. Duesbery and others**' have pointed out that 
stresses which do not result in driving forces on a disloca- 
tion still may affect the overall value of the Peierls stress 
needed to drive the dislocation, because such stresses may 
modify the dislocation core structure,"*^ an effect not ac- 
counted in linear elasticity theory. This effect of non- 
driving stresses is, in fact, one of the common violations 
of the Schmid law which bcc metals exhibit. Diagonal 
stresses, however, do not have a large effect^ on the value 
for the Peierls stress, as they tend to only compress or 
expand the core. Fortunately, because the bulk lattice 
vectors should be relatively close to the that of the dis- 
location cell, we expect all of these effects to be quite 
small, as we have found previously^ and again verify in 
the test calculations below. 

The second simplification in our procedure is that 
rather than applying a stain which imposes a pure Uxz 
stress, we apply a pure e^z stain, which also generates a 
residual axy stress.^" To generate a pure stress of the 
form (Jxz, one would have to apply an additional exy 
stain of a magnitude determined by yet another elastic 
constant of the quadrupolar array. BecaTisc the calcu- 
lation of this constant would increase significantly the 
number of calculations required and because the resid- 
ual ffxy stress^" acts on the plane perpendicular to the 
dislocation, and thus does not create a driving force on 
the dislocations, wc simply apply the pure strain. As 
with the diagonal stress components, although the resid- 
ual in-plane Uxy stress does not drive the dislocations, it 
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can affect the Peierls stress by modifying the core struc- 
ture. Unlike the diagonal stress components, the in-plane 
stress does significantly affect the Peierls stress in bcc 
metals.'*' This effect, however, will be small so long as 
the residual Uxy stress is small compared to the driving 
axz stress. The ratio (7xy|<^xz^ is equal to C" jC where 
C is the elastic constant appearing in Equation (4) and 
C" is another combination of elastic constants. In pure 
bulk cubic materials these constants have the form'^" 

C — -(Cii -f C44 — C12) 
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(2C44 4- C12 — Cii), 



where the Cij are the standard elastic constants for cubic 
materials. 

As evidence of the correlation between the ratio of 
these constants and the errors in the Peierls stress, we 
note that in our previous study^ with same empirical po- 
tential as in Section III, the ratio C" jC varied from 
approximately 1/5 in the smallest cell studied to less 
than 1/10 for all other cells while the error in extract- 
ing the Peierls stress went from 18% to less than 2%, 
respectively. For the potential used in that study, the 
value of C" jC computed from bulk elastic constants is 
1/90. As the bulk value of C" /C within density func- 
tional theory is less than 1/100,^" we expect the errors in 
our ah initio value of the Peierls stress to be even some- 
what smaller. Lending further support to this view is 
the result of Duesbery and Vitek^'^ demonstrating that 
the effect of Gxy stresses on core structure is much less 
for the non-degenerate core structure, which we have in 
our density functional calculations, than for the degen- 
erate core structures, which we had in our inter-atomic 
potential calculations.^ 



B Demonstration 

To demonstrate the efficacy of the above procedure, 
we now proceed to extract the Peierls stress in the (112)- 
direction for vanadium and tantalum from calculations in 
small periodic cells, when the maximum resolved shear 
stress is along a {110}-plane. In all calculations we calcu- 
late the Peierls stress for an infinite straight dislocation 
as our periodic cells has lattice vector 03 = a[lll]/2 along 
the dislocation line. To allow comparison with the Peierls 
stress of isolated dislocations, we employ empirical po- 
tentials for this demonstration. For vanadium we use the 
Finnis-Sinclair''^ potential with modifications made Ack- 
land and Thetford,^^ and for tantalum we use the same 
potential as in Section III but with a slight adjustment 
of parameters to produce a non-degenerate core struc- 
ture. The bulk ratios for C" jC for vanadium and tanta- 
lum within these models are 1/10 and 1/6.5, respectively, 
much larger than the density functional theory value. 

To determine the reference value for the Peierls stress, 
we employed cylindrical boundary conditions with large 



amounts of material, increasing the radius of the cylin- 
ders until the boundary forces were small"'^^ and the 
Peierls stress approached an asymptotic value. To ex- 
tract the Peierls stress from within periodic boundary 
conditions, wc follow precisely the procedure which Sec- 
tion IV A outlines. 
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Strain txz 



FIG. 4: Energy versus strain of vanadium in a quadrupo- 
lar cell of size 42 Ax 20A. Diamonds denote the energy 
calculated at various strains. The line is a quadratic fit 
up to the Peierls stress. 



Figure 4 shows the resulting energy versus stain curve 
for vanadium within a periodic cell of size 42Ax 20A. 
At a strain « 0.054, the curve exhibits a discontinuity 
signaling the critical strain for moving the dislocation. 
The curvature of the fit determines the elastic constant 
C" through Equation (4). Finally, combining this value 
of C with the observed critical strain Equation (5) yields 
the Peierls stress. We repeated this procedure for tanta- 
lum as well. 

Table II summarizes our results for both vanadium 
and tantalum. The table shows that the errors are rela- 
tively small, much smaller than the general discrepancy 
between empirical potentials and the experimental ex- 
trapolations. It is also noteworthy that the potentials 
employed in this demonstration exhibit C" /C ratios over 
an order of magnitude larger than density functional the- 
ory. From these and previous results,^ we conservatively 
estimate that the error in the Peierls stress in the density 
functional calculations below should be no greater than 
« 20%. 





23Ax 11. 5A 


42 Ax 20 A 


Ta 


25% 


10% 


V 


26% 


11% 



TABLE II: Magnitude of percentage error in calculating 
the Peierls stress in periodic cells of two different sizes 
from empirical potentials for vanadium and tantalum. 
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V DENSITY FUNCTIONAL RESULTS AND 
DISCUSSION 

A Computational Details 

All of our first principles electronic structure calcu- 
lations employ the plane-wave density-functional pseu- 
dopotential approach'^'^ within the local density ap- 
proximation.*^''*^ We employ pseudopotential of the 
Kleinman-Bylander form*^ with s, p and d non-local 
channels which has been used successfully in previous 
works^"'"'^ and a plane wave basis with cutoff of 40 Ry- 
dberg. As justified above, we employ a super-cell con- 
taining a quadrupolar array of dislocations of size ri ~ 
5a[l,l,0], r2 = (3/2)a[l,l,2], and = a[l, 1, l]/2, 
where a = 3.25 A is the lattice constant of the cu- 
bic unit cell. The lattice vectors of this cell are ai — 
fi/2 -f2+ r3/2, 02 - ri/2 + f 2 + r3/2 and 03 = ^3. To 
carry out the integrations over the Brillouin zone we use 
a non-zero electronic temperature of ksT = 0.1 eV to 
facilitate integration over the Fermi surface and sample 
the zone at sixteen special fc-points.^*' These choices give 
energy differences reliably to within 0.1 eV/atom. Fi- 
nally, to determine the electronic structure, we minimize 
using the analytically continued functional approach,^* 
expressed within the DFT-|--|-^^ formalism. 



pies electronic structure studies, we compare the energy 
landscape for a dislocation moving along a reaction coor- 
dinate from the easy core configuration to the hard core 
configuration'^' for both molybdenum and tantalum as 
determined within MGPT^"^ and as calculated ab initio. 
Within MGPT, we have carried out the calculation for 
molybdenum ourselves and we used the hard-easy core 
energy difference reported in the literature*^ for tanta- 
lum. Within density functional theory, we have calcu- 
lated the energy at a number of points along the reaction 
pathway for tantalum, and for molybdenum we report the 
difference between the easy core and hard core configu- 
rations as found in Reference 20. We also note that, for 
molybdenum, within density functional theory, the hard 
core configuration was not stable and therefore the stable 
structure found within MGPT was used as the reference 
state. 

Figure 5 shows the results. Most noticeably, the atom- 
istic landscapes are three times stiffer than the ab initio 
landscapes. This raises the question whether the approx- 
imate factor of three overestimate of theoretical calcula- 
tions over the extrapolation of the experimental Peierls 
stresses to zero-temperature is due to defects in the inter- 
atomic potentials or to failures in the connection between 
the experiments and the theoretical calculations. 
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FIG. 5: Energy landscape per Burgers vector along a 
reaction coordinate when going from the easy core to 
the hard core configurations in tantalum and molybde- 
num: density functional theory (squares in upper pan- 
els), MGPT (curves in lower panels). Ab initio results 
for molybdenum are from Reference 20, and the MGPT 
results for tantalum arc from Reference 12. 



B Energy landscapes 

To demonstrate the discrepancies that occurs between 
the predictions of empirical potentials and first princi- 



C Verification of cell size 



To compute the Peierls stress, we shall employ a cell 
of dimensions 23Ax 12A. To verify that long-range elec- 
tronic structure effects in metals do not interfere with 
results in such a cell, we compare the core structure re- 
ported previously^" for this cell with a new calculation 
using a larger cell. Figure 6 shows the result for the core 
structure in cell of size 4lAx 20A. The quadrupole array 
has size ri = 9a[l,i,0], r2 — (5/2)a[l, 1, 2], and rs = 
a[l, 1, l]/2. The lattice vectors are ai ~ fi/2 — r2 + 
7^3/2, 0,2 = n./2 + r2 + Tzl'i. and 03 = fa. Here we use 
eight special A:-points,^'^ which is sufficient as the lattice 
vectors in the plane of the dislocation have doubled. 

We note that the core structure is very similar to pre- 
vious results^*^ which used a smaller cell equal in size to 
that we shall employ for our calculation of the Peierls 
stress. We therefore do not expect the long-range nature 
of electronic effects in metallic systems to greatly affect 
the value which shall extract for the Peierls stress. We 
also note that the empirical potential results for tantalum 
(Section IV) had a very large cutoff of 9A and accurate 
results were obtained in the smallest cells used in those 
calculations. These facts lend confidence in the reliability 
of our density-functional theory predictions below. 
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FIG. 6: Easy core structure for tantalum calculated 
within density functional theory. This cell of size 41 Ax 
20 A gives very similar results to the cell of size 23A X 
12 A, used in Reference 20. The solid triangle in the figure 
represents the center of the dislocation. 
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FIG. 8: Differential displacement map at an applied 
strain of 0.04. {axz = 1-41 GPa) The cell has converged 
to within the tolerance defined in the text. The solid 
triangle in the figure is the location of the center of the 
dislocation, under no stress. The center of the disloca- 
tion, at this applied stress, has not moved. 



D Results for the Peierls Stress 

Figure 7 shows our density functional theory results for 
energy as a function of strain following the procedure of 
Section IV A. We regard each data point as fully relaxed 
when the magnitude of the residual force on each atom 
is less than 0.005 eV/A. 
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FIG. 7: Energy versus strain for tantalum in a cell of size 
23Ax 12 A, calculated within density functional theory. 
Squares denote the calculate energy at various strains. 
The line is a quadratic fit to the first five data points. 
The graph indicates that the points with strain at and 
above 0.05 {axz = 1-76 GPa) will be at a stress above the 
Peierls stress. 
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FIG. 9: Differential displacement map at an applied 
strains of 0.05. {axz = 1.76 Gpa) The cell has converged 
to within the tolerance defined in the text. This fig- 
ure shows that the core has moved from the old triad 
(solid triangle) to a new triad (upside-down, solid trian- 
gle) along a (112)-direction. 



From Figure 7 it is apparent that, at a strain of 0.05 
the Peierls stress has been exceeded, whereas at a strain 
of 0.04 the energy curve lies on the elastic solution. From 
the curvature of the data in the elastic region (C = 35.2 
GPa) and these critical strains, we bound the Peierls 
stress {axz) to be between 1.41 GPa and 1.76 GPa. These 
results are slightly lower than the previous density func- 
tional results of « 1.8 GPa obtained using Greens func- 
tion boundary conditions.® However both results are 
still in reasonable agreement, as we expect our results 
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to be within 20% of the infinite ceU limit, while it is un- 
clear how the domain boundary affected the value of the 
Peierls stress in the Greens function boundary condition 
calculation. 

Figures 8 and 9 show differential displacement maps"^ 
of the dislocation configuration at strains of 0.04 and 
0.05, respectively. The solid triangles in the figures rep- 
resent the location of the center of the dislocation before 
the application of stain and the upside-down, solid tri- 
angle represents the center of the dislocation once it has 
moved. Figure 8 shows the relaxed dislocation core to 
remain in the position of its unstressed state, whereas 
Figure 9 shows that at an applied strain of 0.05, the dis- 
location moves onto the next triad of columns of atoms. 
This glide is consistent with such screw dislocations as 
the initial displacement is along the (Oll)-plane. The 
subsequent motion should be along the (liO)-plane, and 
hence overall motion is along a {112}-plane (twinning di- 
rection), consistent with the previous density functional 
theory^ calculations. 

Intriguingly, despite the fact that the ab initio energy 
landscape is less corrugated than that of the inter-atomic 
potential by a factor of nearly three (Figure 5) , the above 
ab initio result for the Peierls stress is over a factor of 
two larger than the empirical potential result (Table I). 
Moreover, our result is over a factor of five larger than ex- 
trapolations of experimental data to zero-temperature, "'^^ 
a discrepancy much larger than any effects from the use 
of either our relatively small quadrupolar cell or the lo- 
cal density approximation to density-functional theory. 
Certainly, as Section III shows, some of this discrep- 
ancy can be due to the extrapolation of experimental 
data to zero temperature. However, errors in the zero- 
temperature extrapolation may not likely account for all 
of the discrepancy as the underestimation in Figure 2 is 
only Ri 10 — 20%. One must therefore consider the possi- 
bility of other factors such as the effects of kinks, defects, 
surfaces, strain-rate dependencies and screening from the 
presence of other dislocations, to accurately predict the 
experimental results. 



VI CONCLUSION 

This work explores various aspects of the the Peierls 
stress for the (111) screw dislocation in bcc tantalum. 
The first non-zero temperature results for the Peierls 
stress in this system shows both a strong orientation and 
temperature dependent response, consistent with experi- 
mental results. This data also demonstrate that common 
extrapolations of experimental data tend to underesti- 
mate the zero-temperature limit. 

We have also presented the first density functional 
theory calculation for the Peierls stress within perioidic 
boundary conditions, the approach best suited to metal- 
lic systems. The value we find for the Peierls stress is 
substantially larger than both the experimental extrapo- 



lations and current empirical potential results. This dif- 
ference is much larger than errors which could normally 
be attributed to the use of a relatively small unit cell 
or the local density approximation to density-functional 
theory. The error is also significantly larger than the 
underestimation we have seen in the the extrapolation 
of nonzero-temperature data to zero temperature, thus 
supporting the notion that mechanisms other than the 
simple Peierls resistance may play an important role in 
controlling the process of low temperature slip. 
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